

# Hicks, William D. (hickswd@appstate.edu). 2015. "Partisan Competition 
# and the Efficiency of Lawmaking in American Legislatures, 1991-2009"
# American Politics Research 43(5): 743-770. 


# This code, if run properly, replicates Table 1, and Figurres 1-3 #
# Please cite the above article if you use these data or code #
 
 
# Required packages

library(nlme)
library(Hmisc)

#import .csv 
efficiency <- read.csv("~REPLICATION_DATA.csv")# write your own path here



#VARIABLE LABELS AND DEFINITIONS 
label(efficiency$statecode)<-"Unique code for each state"
label(efficiency$year)<-"Year"
label(efficiency$enacbills)<-"Number of enacted bills in legislature i, year t"
label(efficiency$intbills)<-"Number of introduced bills in legislature i, year t"
efficiency$intbills00<-efficiency$intbills/100
label(efficiency$intbills00)<-"Hundreds of bill introductions in legislature i, year t"
label(efficiency$leg.comp)<-"Folded parstisan seat margin in legislature i, year t"
label(efficiency$polar)<-"Polarization in legislature i, year t"
label(efficiency$divided_gov)<-"Divided government in state i, year t"
label(efficiency$term.limits.yoi)<-"Term Limits (adopted & implemented) in legislature i, year t"
label(efficiency$profess.bbs)<-"Legislative professionalism in legislature i, year t"
label(efficiency$compbbs_emaint)<-"Mean margin of victory for legislators in legislature i, year t"
label(efficiency$bicam.diffs.floor)<-"Bicameral differences (i.e., ideological distance between median House & Senate members) in legislature i, year t"
label(efficiency$time)<-"Time (year=time+1991)"
label(efficiency$group.index)<-"Interest group power index for state i, year t"
label(efficiency$pop.size.mil)<-"Population in millions for state i, year t"
label(efficiency$racial.index)<-"Racial diversity index for state i, year t"



#CENTERING VARIABLES FOR MODELS

efficiency$cmean.legcomp<-efficiency$leg.comp-mean(efficiency$leg.comp, na.rm=TRUE)
efficiency$cmean.polar<-efficiency$polar-mean(efficiency$polar, na.rm=TRUE)
efficiency$cmean.legcompXpolar<-efficiency$cmean.polar*efficiency$cmean.legcomp
efficiency$cmean.legcompXdivgov<-efficiency$cmean.legcomp*efficiency$divided_gov
efficiency$cmean.bicam.diffs.floor<-efficiency$bicam.diffs.floor-mean(efficiency$bicam.diffs.floor, na.rm=TRUE)
efficiency$cmean.compbbs_emaint<-efficiency$compbbs_emaint-mean(efficiency$compbbs_emaint, na.rm=TRUE)
efficiency$cmean.profess.bbs<-efficiency$profess.bbs-mean(efficiency$profess.bbs, na.rm=TRUE)
efficiency$cmean.intbills00<-efficiency$intbills00-mean(efficiency$intbills00, na.rm=TRUE)
efficiency$cmean.group.index<-efficiency$group.index-mean(efficiency$group.index, na.rm=TRUE)
efficiency$cmean.pop.size.mil<-efficiency$pop.size.mil-mean(efficiency$pop.size.mil, na.rm=TRUE)
efficiency$cmean.racial.index<-efficiency$racial.index-mean(efficiency$racial.index, na.rm=TRUE)


#################################################################################### 
#################################################################################### 
#APR TABLE 1 
####################################################################################
####################################################################################


M.FULL<-lme(enacbills ~ cmean.legcomp + cmean.polar + divided_gov + cmean.bicam.diffs.floor + cmean.compbbs_emaint + cmean.profess.bbs + term.limits.yoi + divided_gov + time + cmean.intbills00 + I(cmean.intbills00^2) + cmean.legcompXpolar + cmean.legcompXdivgov + cmean.group.index + cmean.pop.size.mil + cmean.racial.index , random=~time|state,correlation=corAR1(), method="ML",  data=efficiency, na.action=na.omit)

round(summary(M.FULL)$tTable, dig=4)
AIC(M.FULL)
BIC(M.FULL)
anova(M.FULL)
VarCorr(M.FULL)


#######################################################################################
#FUNCTIONS TO CONSTRUCT MARGINAL EFFECTS FOR INTERACTIONS AND POLYNOMIALS (FIGURES 2 & 3)
#BOTH REQUIRE THAT I USE THE MODELS WITH COVARIATES IN A PARTICULAR ORDER
#FOR ME.f(M, X, Z, xlab, ylab, level), lme(y~x+z+c1+...+ck+xz)
#FOR POLY.ME(M, X, xlab, ylab, level), lme(y~x+x2+c1+...ck)  
#######################################################################################



ME.f<-function(M,X,Z,xlab,ylab,level)
{
S <- summary(M)
N <- c(1:20)

#     ################################################################  # #       Create 20 values in a vector between min                        # #       and max on the Z variable                                       # #     ################################################################  #

zmin <- rep(min(Z,na.rm=TRUE), 20)
zmax <- rep(max(Z, na.rm=TRUE), 20)
znew <- (((N-1)/(20-1))*(zmax-zmin))+zmin

#     ################################################################  # #       Grab elements of the coefficient and variance-covariance matrix # #       to calculate the marginal effect and standard                   # #       confidence intervals.                                           # #     ################################################################  #


H <- head(fixef(S),3)
T <- tail(fixef(S),1)
b <- rbind(H,T)
b<-as.data.frame(b)

Vcov <- vcov(M)
Vcov <- as.data.frame(Vcov)
Vcov1 <- Vcov[,c(1:3)]
Vcov2 <- Vcov[,-c(0:0-length(Vcov))]
Vcov <- cbind(Vcov1,Vcov2)

Vh <- head(Vcov,3)
Vt <- tail(Vcov,1)
V <- rbind(Vh,Vt)

b1 <- b[1,2]
b2 <- b[1,3]
b3 <- b[2,1]

varb1 <- V[2,2]
varb2 <- V[3,3]
varb3 <- V[4,4]

covb1b3 <- V[4,2]
covb2b3 <- V[4,3]

#     ################################################################ # #      The following generates the marginal effect of x on y           #
 #      for values of z                                                 # #     ################################################################ #

conb <- b1+b3*znew

#     ################################################################  # #       Following calculates the s.e.'s for the marginal effect of X    # #       for all values of Z                                             # #     ################################################################  #

conse <- sqrt(varb1 + varb3*((znew)^2) + 2*covb1b3*znew)

#     ################################################################  # #       Upper and Lower confidence intervals using                     #
#        approximations of 1.95 for 95 and 1.645 for 90                 #          #     ################################################################  #

ci<-NA
ci[level==95]<-1.96
ci[level==90]<-1.645

a=ci*conse
upper=conb+a
lower=conb-a

#     ################################################################  # #       Graph the marginal effect of X on Y across the desired range of # #       the modifying variable Z.                                       # #     ################################################################  #
plot(c(znew,znew,znew), c(a,upper,lower), type="n",xlab=xlab,ylab=ylab)+rug(jitter(Z))

polygon(c(znew,rev(znew)),c(lower,rev(upper)),col = "grey85", border =NA)
lines(znew,conb,type="l", lwd=2)
abline(h=0, lty=1, lwd=0.5)
}

#############REPLACE POLYGON-LINES TO USE DIFFERENT CONFIDENCE INTERVALS###############
#lines(znew,conb,col="black")
#lines(znew,upper,col="black",lty=2)
#lines(znew,lower,col="black",lty=2)
######################################################################################

######################################################################################
######################################################################################


POLY.ME<-function(M,X,xlab,ylab,level)
{
S <- summary(M)
N <- c(1:20)

#     ################################################################  # #       Create 20 equally spaced values in a vector between min         # #       and max on the X variable                                       # #     ################################################################  #

xmin <- rep(min(X,na.rm=TRUE), 20)
xmax <- rep(max(X, na.rm=TRUE), 20)
xnew <- (((N-1)/(20-1))*(xmax-xmin))+xmin

#     ################################################################  # #       Grab elements of the coefficient and variance-covariance matrix # #       to calculate the marginal effect and standard                   # #       confidence intervals.                                           # #     ################################################################  #


H <- head(fixef(S),3)
b<-as.data.frame(H)

Vcov <- vcov(M)
Vcov <- as.data.frame(Vcov)

Vh <- head(Vcov,3)
Vt <- tail(Vcov,1)
V <- rbind(Vh,Vt)

b1 <- b[2,1]
b2 <- b[3,1]

varb1 <- V[2,2]
varb2 <- V[3,3]
covb1b2 <- V[3,2]


#     ################################################################  # #       Calculate the marginal effect of X on Y for all                 # #       values of *X*.                                                   # #     ################################################################  #

conb <- b1+(2*b2*xnew)

#     ################################################################  # #       Calculate the standard errors for the marginal effect of X on Y # #     ################################################################  #

conse <- sqrt(varb1 + 4*varb2*((xnew)^2) + 4*covb1b2*xnew)

#     ################################################################  # #              Generate Upper and Lower confidence intervals            #         
#              using  1.96 for 95 1.645 for 90                           # #     ################################################################  #

ci<-NA
ci[level==95]<-1.96
ci[level==90]<-1.645

a=ci*conse
upper=conb+a
lower=conb-a

#     ################################################################  # #       Graph the marginal effect of X on Y across the desired range of # #                                x                                       # #     ################################################################  #
plot(c(xnew,xnew,xnew), c(a,upper,lower), type="n",xlab=xlab,ylab=ylab)+rug(jitter(X))

polygon(c(xnew,rev(xnew)),c(lower,rev(upper)),col = "grey85", border =NA)
lines(xnew,conb,type="l", lwd=2)
abline(h=0, lty=1, lwd=0.5)
}

##################################################################################
################################################################################## 

#################################################################################### 
#################################################################################### 
#APR FIGURE 1 (I DO THIS ONE WITHOUT THE FUNCTIONS, OTHERS WITH)
####################################################################################
####################################################################################

M.FULL.FIGURE<-lme(enacbills ~ cmean.legcomp + cmean.polar + divided_gov + cmean.bicam.diffs.floor + cmean.compbbs_emaint + cmean.profess.bbs + term.limits.yoi + divided_gov + time + cmean.intbills00 + I(cmean.intbills00^2) + cmean.group.index + cmean.pop.size.mil + cmean.racial.index + cmean.legcompXpolar + cmean.legcompXdivgov, random=~time|state,correlation=corAR1(), method="ML",  data=efficiency, na.action=na.omit)

S <- summary(M.FULL.FIGURE)
N <- c(1:20)
zmin <- rep(min(efficiency$cmean.polar,na.rm=TRUE), 20)
zmax <- rep(max(efficiency$cmean.polar, na.rm=TRUE), 20)
znew <- (((N-1)/(20-1))*(zmax-zmin))+zmin
wnew1<-min(efficiency$divided_gov, na.rm=TRUE)
wnew3<-max(efficiency$divided_gov, na.rm=TRUE)
H <- head(fixef(S),4)
T <- tail(fixef(S),2)
b <- rbind(H,T)
b<-as.data.frame(b)
Vcov <- vcov(M.FULL.FIGURE)
Vcov <- as.data.frame(Vcov)
Vcov1 <- Vcov[,c(1:4)]
Vcov2 <- Vcov[,-c(1:0-length(Vcov))]
Vcov <- cbind(Vcov1,Vcov2)
Vh <- head(Vcov,4)
Vt <- tail(Vcov,2)
V <- rbind(Vh,Vt)
b1 <- b[1,2]
b2 <- b[1,3]
b3 <- b[1,4]
b4 <- b[2,1]
b5 <- b[2,2]
varb1 <- V[2,2]
varb2 <- V[3,3]
varb3 <- V[4,4]
varb4 <- V[5,5]
varb5 <- V[6,6]
covb1b4 <- V[2,5]
covb1b5 <- V[2,6]
covb4b5 <- V[5,6]
conb1 <- b1+b4*znew+b5*wnew1
conse1 <- sqrt(varb1 + varb4*((znew)^2)+varb5*((wnew1)^2)+(2*covb1b4*znew)+(2*covb1b5*wnew1)+(2*znew*wnew1*covb4b5))
conb3 <- b1+b4*znew+b5*wnew3
conse3 <- sqrt(varb1 + varb4*((znew)^2)+varb5*((wnew3)^2)+(2*covb1b4*znew)+(2*covb1b5*wnew3)+(2*znew*wnew3*covb4b5))
ci<-1.645
a1=ci*conse1
a3=ci*conse3
upper1=conb1+a1
lower1=conb1-a1
upper3=conb3+a3
lower3=conb3-a3

par(mfrow=c(1,2))
plot(znew, conb1, xlab="Polarization", ylab="ME(Folded Partisan Seat Margin | Unified Gov)", type="n", xlim=c(min(efficiency$cmean.polar, na.rm=T),max(efficiency$cmean.polar, na.rm=T)), ylim=c(-15,10))
rug(jitter(efficiency$cmean.polar))
polygon(c(znew,rev(znew)),c(lower1,rev(upper1)),col = "grey85", border=NA)
lines(znew,conb1,type="l", lwd=2)
abline(h=0, lty=1, lwd=0.5)

plot(znew, conb3, xlab="Polarization", ylab="ME(Folded Partisan Seat Margin | Divided Gov)", type="n", xlim=c(min(efficiency$cmean.polar, na.rm=T),max(efficiency$cmean.polar, na.rm=T)), ylim=c(-15,10))
rug(jitter(efficiency$cmean.polar))
polygon(c(znew,rev(znew)),c(lower3,rev(upper3)),col = "grey85", border=NA)
lines(znew,conb3,type="l", lwd=2)
abline(h=0, lty=1, lwd=0.5)

#################################################################################### 
####################################################################################
#APR FIGURE 2
#################################################################################### 
#################################################################################### 


M.FULL.FIGURE21<-lme(enacbills ~ cmean.polar + cmean.legcomp + divided_gov + cmean.bicam.diffs.floor + cmean.compbbs_emaint + cmean.profess.bbs + term.limits.yoi + divided_gov + time + cmean.intbills00 + I(cmean.intbills00^2) + cmean.group.index + cmean.pop.size.mil + cmean.racial.index + cmean.legcompXdivgov+ cmean.legcompXpolar, random=~time|state,correlation=corAR1(), method="ML",  data=efficiency, na.action=na.omit)
M.FULL.FIGURE22<-lme(enacbills ~ divided_gov + cmean.polar + cmean.legcomp +  cmean.bicam.diffs.floor + cmean.compbbs_emaint + cmean.profess.bbs + term.limits.yoi + divided_gov + time + cmean.intbills00 + I(cmean.intbills00^2) + cmean.group.index + cmean.pop.size.mil + cmean.racial.index + cmean.legcompXpolar + cmean.legcompXdivgov, random=~time|state,correlation=corAR1(), method="ML",  data=efficiency, na.action=na.omit)
par(mfrow=c(1,2))
ME.f(M.FULL.FIGURE21, efficiency$cmean.polar,efficiency$cmean.legcomp,"Folded Partisan Seat Margin","Marginal Effect of Polarization",90)
ME.f(M.FULL.FIGURE22, efficiency$divided_gov,efficiency$cmean.legcomp,"Folded Partisan Seat Margin","Marginal Effect of Divided Government",90)


#################################################################################### 
####################################################################################
#APR FIGURE 3
#################################################################################### 
#################################################################################### 

M.FULL31<-lme(enacbills ~ cmean.intbills00 + I(cmean.intbills00^2) + cmean.legcomp + cmean.polar + divided_gov + cmean.bicam.diffs.floor + cmean.compbbs_emaint + cmean.profess.bbs + term.limits.yoi + divided_gov + time + cmean.legcompXpolar + cmean.legcompXdivgov + cmean.group.index + cmean.pop.size.mil + cmean.racial.index , random=~time|state,correlation=corAR1(), method="ML",  data=efficiency, na.action=na.omit)

#MEAN CMEAN.BILL INTRO IN HUNDREDS IS 0, MIN IS -19.40, MAX IS 162

par(mfrow=c(1,2))
quadfit<-function(x) 503.2305+7.8020*x-0.0359*x^2
plot(quadfit, -19, 160, ylab="Expected Number of Bills Enacted", xlab="Hundreds of Bill Introductions")
POLY.ME(M.FULL31, efficiency$cmean.intbills00, xlab="Hundreds of Bill Introductions", ylab="Marginal Effect of Bill Introductions", 90)




